function res = vmr_InhomogeneityCheck(hfile, tal, acpc)
% VMR::InhomogeneityCheck  - perform inhomogeneity check on VMR
%
% FORMAT:       [res = ] vmr.InhomogeneityCheck([tal [, acpc]])
%
% Input fields:
%
%       tal         optional TAL information, if VMR in ACPC/native space
%                   (filename or object)
%       trf         optional ACPC transformation, if VMR in native space
%                   (filename or object)
%
% Output fields:
%
%       res         1x1 struct with fields
%        .meanx/y/z mean intensities of by prior selected voxels
%        .numx/y/z  number of selected voxels
%        .voi       VOI struct with coordinates of selected voxels (TAL)
%                   suitable to add to a VOI object
%
% Note: this function uses a prior file compiled out
%       20 WM segmentations (0.5mm coordinate list)
%
% Note2: the mean/numx/y/z are in BV's internal logic !

% Version:  v0.8a
% Build:    9110918
% Date:     Nov-09 2009, 6:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'vmr')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if bc.FramingCube > 256
    error( ...
        'BVQXfile:BadArgument', ...
        'Only valid for up to 256 FramingCube.' ...
    );
end
lobj = cell(1, 3);
if nargin < 2 || ...
   ((numel(tal) ~= 1 || ...
     ~isBVQXfile(tal, 'tal')) && ...
    (~ischar(tal) || ...
     numel(tal) < 5 || ...
     ~strcmpi(tal(end-3:end), '.tal') || ...
     exist(tal(:)', 'file') ~= 2))
    tal = {};
elseif ischar(tal)
    try
        tal = BVQXfile(tal(:)');
        lobj{1} = tal;
    catch
        rethrow(lasterror);
    end
else
    tal = {tal};
end
if nargin < 3 || ...
   ((numel(acpc) ~= 1 || ...
     ~isBVQXfile(acpc, 'trf')) && ...
    (~ischar(acpc) || ...
     numel(acpc) < 5 || ...
     ~strcmpi(acpc(end-3:end), '.trf') || ...
     isempty(regexpi(acpc(:)', '^.*acpc.*$')) || ...
     exist(acpc(:)', 'file') ~= 2))
    acpc = {};
elseif ischar(acpc)
    try
        acpc = BVQXfile(acpc(:)');
        lobj{2} = acpc;
    catch
        clearbvqxobjects(lobj);
        rethrow(lasterror);
    end
else
    acpc = {acpc};
end

% load coordinate list
try
    wm_cs = load([bvqxtools_path('ihc') '/prob_20_tal_wm_cs.mat']);
    wm_cs = wm_cs.wm_cs;
catch
    clearbvqxobjects(lobj);
    error( ...
        'BVQXfile:InternalError', ...
        'Error loading prior probability coordinates.' ...
    );
end

% create temporary VOI
voi = BVQXfile('new:voi');
lobj{3} = voi;
voic = bvqxfile_getcont(voi.L);
voic.VOI(1).Name = 'WM';
voic.VOI = voic.VOI(1);
voic.VOI.Color = [255 255 255];
voic.VOI.Voxels = double(wm_cs);
voic.VOI.NrOfVoxels = size(wm_cs, 1);
if ~isempty(tal) || ...
    ~isempty(acpc)
    bvqxfile_setcont(voi.L, voic);
    voi_ApplyTrf(voi, {tal{:}, acpc{:}}, true, true);
    voic = bvqxfile_getcont(voi.L);
end

% prepare output
res = struct;
res.meanx = NaN(1, 256);
res.meany = NaN(1, 256);
res.meanz = NaN(1, 256);
res.medianx = NaN(1, 256);
res.mediany = NaN(1, 256);
res.medianz = NaN(1, 256);
res.numx  = zeros(1, 256);
res.numy  = zeros(1, 256);
res.numz  = zeros(1, 256);
res.voi = voic.VOI;

% clear VOI object
bvqxfile_clear(voi.L);
lobj{3} = [];

% get BV coordinates
bvc = 128 - unique(round(double(voic.VOI.Voxels(:, [2, 3, 1]))), 'rows');

% deduce offset
bvc(:, 1) = 1 + bvc(:, 1) - (bc.OffsetX + 0.5 * (256 - ceil(bc.FramingCube)));
bvc(:, 2) = 1 + bvc(:, 2) - (bc.OffsetY + 0.5 * (256 - ceil(bc.FramingCube)));
bvc(:, 3) = 1 + bvc(:, 3) - (bc.OffsetZ + 0.5 * (256 - ceil(bc.FramingCube)));

% what data
if isequal(size(bc.VMRData), size(bc.VMRData16))
    vd = bc.VMRData16(:, :, :);
else
    vd = bc.VMRData(:, :, :);
end

% iterate over axes
vs = size(vd);
bvc((bvc(:, 1) < 1 | bvc(:, 1) > vs(1)), :) = [];
bvc((bvc(:, 2) < 1 | bvc(:, 2) > vs(2)), :) = [];
bvc((bvc(:, 3) < 1 | bvc(:, 3) > vs(3)), :) = [];
bvc = sub2ind(vs, bvc(:, 1), bvc(:, 2), bvc(:, 3));
vds = vd(bvc);
vd(:) = 0;
vd(bvc) = vds;
bvc = (vd > 0);
for c = 1:vs(1)
    vds = vd(c, :, :);
    vds = vds(bvc(c, :, :));
    res.numx(c) = numel(vds);
    if res.numx(c) > 1
        res.meanx(c) = mean(vds);
        res.medianx(c) = median(vds);
    end
end
for c = 1:vs(2)
    vds = vd(:, c, :);
    vds = vds(bvc(:, c, :));
    res.numy(c) = numel(vds);
    if res.numy(c) > 1
        res.meany(c) = mean(vds);
        res.mediany(c) = median(vds);
    end
end
for c = 1:vs(3)
    vds = vd(:, :, c);
    vds = vds(bvc(:, :, c));
    res.numz(c) = numel(vds);
    if res.numz(c) > 1
        res.meanz(c) = mean(vds);
        res.medianz(c) = median(vds);
    end
end

% clear objects
clearbvqxobjects(lobj);
